miχpods: MOM6#

  • baseline KPP missing

  • bootstrap error on mean, median?

  • Try daily vs hourly

  • add TAO χpods

  • fix EUC max at 110

  • figure out time base

  • check eps vs Ri histogram

    • I think my Ri is too low?

Questions:

  • N2 vs N2T

  • match time intervals

  • match vertical resolution and extent

  • restrict TAO S2 to ADCP depth range

  • restrict to top 200m.

Notes:

  • Filtering out MLD makes a big difference!

  • SInce we’re working with derivatives does the vertical grid matter (as long as it is coarser than observations)?

1 difvho ocean_vertical_heat_diffusivity 
2 difvso ocean_vertical_salt_diffusivity

  • Need lateral / neutral terms for total χ, ε

  • Why can’t I reproduce the ε figure?

  • Why can’t I save the new station data

  • Need to consider Ri smoothing in models: From Gustavo:

    You mentioned some kind of inconsistency between your diagnostics and what MOM6 was doing for the interior shear mixing scheme. I am working to implement the option to apply vertical smoothing in Ri multiple times, instead of only once which is what we are doing right now. I noticed that the diagnostic ri_grad_shear is saved before the smoothing is applied. There is another diagnostic (ri_grad_shear_smooth) that saves the smoothed Ri. Perhaps you were looking at ri_grad_shear instead of ri_grad_shear_smooth and this can explain the inconsistency.

Diagnostics notes#

  • go through prepare; pass in dataset and optional xgcm.Grid

  • different pre-processing steps for different datasets

  • order, sign of z coordinate is painful; here normalizing

  • composing with matplotlib subfigures

References#

Warner & Moum (2019)

Setup#

%load_ext watermark

import datetime
import glob
import os
import warnings

import cf_xarray as cfxr
import dask
import dask_jobqueue
import dcpy
import distributed
import flox.xarray
import holoviews as hv
import hvplot.xarray
import matplotlib as mpl
import matplotlib.pyplot as plt
import mom6_tools
import numpy as np
import pandas as pd
import xarray as xr
import xgcm
from holoviews import opts

%aimport pump
from pump import mixpods

hv.notebook_extension("bokeh")

dask.config.set({"array.slicing.split_large_chunks": False})
plt.style.use("bmh")
plt.rcParams["figure.dpi"] = 140
xr.set_options(keep_attrs=True, display_expand_data=False)

gcmdir = "/glade/campaign/cgd/oce/people/bachman/TPOS_1_20_20_year/OUTPUT/"  # MITgcm output directory
stationdirname = gcmdir

mixing_layer_depth_criteria = {
    "boundary_layer_depth": {"name": "KPPhbl|KPP_OBLdepth|ePBL_h_ML"},
}

%watermark -iv

/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask_jobqueue/core.py:20: FutureWarning: tmpfile is deprecated and will be removed in a future release. Please use dask.utils.tmpfile instead.
  from distributed.utils import tmpfile

dcpy         : 0.1.dev385+g121534c
dask         : 2022.11.1
sys          : 3.10.8 | packaged by conda-forge | (main, Nov 22 2022, 08:26:04) [GCC 10.4.0]
dask_jobqueue: 0.7.3
pump         : 0.1
pandas       : 1.5.2
json         : 2.0.9
flox         : 0.6.4
mom6_tools   : 0.0.post244
xarray       : 2022.11.0
xgcm         : 0.6.1
numpy        : 1.23.5
matplotlib   : 3.6.2
holoviews    : 1.15.2
cf_xarray    : 0.7.5
hvplot       : 0.8.2
distributed  : 2022.11.1

if "client" in locals():
    client.close()
    del client
if "cluster" in locals():
    cluster.close()

# env = {"OMP_NUM_THREADS": "3", "NUMBA_NUM_THREADS": "3"}

# cluster = distributed.LocalCluster(
#    n_workers=8,
#    threads_per_worker=1,
#    env=env
# )

if "cluster" in locals():
    del cluster

# cluster = ncar_jobqueue.NCARCluster(
#    project="NCGD0011",
#    scheduler_options=dict(dashboard_address=":9797"),
# )
# cluster = dask_jobqueue.PBSCluster(
#    cores=9, processes=9, memory="108GB", walltime="02:00:00", project="NCGD0043",
#    env_extra=env,
# )


cluster = dask_jobqueue.PBSCluster(
    cores=1,  # The number of cores you want
    memory="23GB",  # Amount of memory
    processes=1,  # How many processes
    queue="casper",  # The type of queue to utilize (/glade/u/apps/dav/opt/usr/bin/execcasper)
    # log_directory="/glade/scratch/dcherian/dask/",  # Use your local directory
    resource_spec="select=1:ncpus=1:mem=23GB",  # Specify resources
    project="ncgd0011",  # Input your project ID here
    walltime="02:00:00",  # Amount of wall time
    interface="ib0",  # Interface to use
)
cluster.adapt(maximum_jobs=24, minimum_jobs=3)
client = distributed.Client(cluster)
client

Client

Client-f9c48d0b-7199-11ed-84ca-3cecef1ac722

Connection method: Cluster object Cluster type: dask_jobqueue.PBSCluster
Dashboard: https://jupyterhub.hpc.ucar.edu/stable/user/dcherian/proxy/8787/status

Cluster Info

Read data#

TAO#

tao_Ri = xr.load_dataarray(
    "tao-hourly-Ri-seasonal-percentiles.nc"
).cf.guess_coord_axis()

tao_gridded = xr.open_dataset(
    os.path.expanduser("~/work/pump/zarrs/tao-gridded-ancillary.zarr"),
    chunks="auto",
    engine="zarr",
).sel(longitude=-140, time=slice("1996", None))
tao_gridded["depth"].attrs["axis"] = "Z"
tao_gridded["shallowest"].attrs.clear()

chipod = (
    dcpy.oceans.read_cchdo_chipod_file(
        "~/datasets/microstructure/osu/chipods_0_140W.nc"
    )
    .sel(time=slice("2015"))
    # move from time on the half hour to on the hour
    .coarsen(time=2, boundary="trim")
    .mean()
    # "Only χpods between 29 and 69 m are used in this analysis as
    # deeper χpods are more strongly influenced by the variability of zEUC than by surface forcing."
    # - Warner and Moum (2019)
    .sel(depth=slice(29, 69))
    .reindex(time=tao_gridded.time, method="nearest", tolerance="5min")
    .pipe(mixpods.normalize_z, sort=True)
)

tao_gridded = tao_gridded.update(
    chipod[["chi", "KT", "eps", "Jq"]]
    .reset_coords(drop=True)
    .rename({"depth": "depthchi"})
)

tao_gridded = mixpods.prepare(tao_gridded, oni=pump.obs.process_oni())

tao_gridded = tao_gridded.update(
    mixpods.pdf_N2S2(tao_gridded.drop_vars(["shallowest", "zeuc"]))
)
tao_gridded = mixpods.load(tao_gridded)

/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/core/dataset.py:247: UserWarning: The specified Dask chunks separate the stored chunks along dimension "depth" starting at index 42. This could degrade performance. Instead, consider rechunking after loading.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/core/dataset.py:247: UserWarning: The specified Dask chunks separate the stored chunks along dimension "time" starting at index 199726. This could degrade performance. Instead, consider rechunking after loading.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/core/dataset.py:247: UserWarning: The specified Dask chunks separate the stored chunks along dimension "longitude" starting at index 2. This could degrade performance. Instead, consider rechunking after loading.
  warnings.warn(

np.log10(tao_gridded.eps).sel(time=slice("2005", "2015")).resample(
    time="M"
).mean().hvplot.quadmesh(clim=(-7.5, -6))

sub = tao_gridded.sel(time="2010-01")

t = sub.theta.hvplot.quadmesh(cmap="turbo_r")
dt = (
    sub.theta - sub.theta.reset_coords(drop=True).cf.sel(Z=[0, -5]).cf.max("Z")
).hvplot.quadmesh(clim=(-0.15, 0.15), cmap="RdBu_r")
newmld = mixpods.get_mld_tao_theta(sub.theta.reset_coords(drop=True))

(
    dt
    * sub.reset_coords().mldT.hvplot.line(color="w", line_width=2)
    * newmld.reset_coords(drop=True).hvplot.line(color="orange", line_width=1)
).opts(frame_width=1200)

(
    tao_gridded.reset_coords().mldT.resample(time="5D").mean().hvplot.line()
    * mixpods.get_mld_tao_theta((tao_gridded.reset_coords().theta))
    .resample(time="5D")
    .mean()
    .hvplot.line()
)

tao_gridded.u.cf.plot()
tao_gridded.eucmax.plot()

[<matplotlib.lines.Line2D at 0x2b38c6a4a8c0>]
_images/mixpods-mom6_11_1.png

MITgcm stations#

stations = pump.model.read_stations_20(stationdirname)

gcmeq = stations.sel(
    longitude=[-155.025, -140.025, -125.025, -110.025], method="nearest"
)
# enso = pump.obs.make_enso_mask()
# mitgcm["enso"] = enso.reindex(time=mitgcm.time.data, method="nearest")

# gcmeq["eucmax"] = pump.calc.get_euc_max(gcmeq.u)
# pump.calc.calc_reduced_shear(gcmeq)

# oni = pump.obs.process_oni()
# gcmeq["enso_transition"] = mixpods.make_enso_transition_mask(oni).reindex(time=gcmeq.time.data, method="nearest")


mitgcm = gcmeq.sel(longitude=-140.025, method="nearest")
metrics = mixpods.normalize_z(pump.model.read_metrics(stationdirname), sort=True)
mitgcm = mixpods.normalize_z(mitgcm, sort=True)

mitgcm_grid = xgcm.Grid(
    metrics.sel(latitude=mitgcm.latitude, longitude=mitgcm.longitude, method="nearest"),
    coords=({"Z": {"center": "depth", "outer": "RF"}, "Y": {"center": "latitude"}}),
    metrics={"Z": ("drF", "drC")},
    periodic=False,
    boundary="fill",
    fill_value=np.nan,
)

mitgcm.theta.attrs["standard_name"] = "sea_water_potential_temperature"
mitgcm.salt.attrs["standard_name"] = "sea_water_salinity"
mitgcm["KPPviscAz"].attrs["standard_name"] = "ocean_vertical_viscosity"

mitgcm = (
    mixpods.prepare(mitgcm, grid=mitgcm_grid, oni=pump.obs.process_oni())
    .sel(latitude=0, method="nearest")
    .cf.sel(Z=slice(-250, 0))
)

/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask/array/core.py:4806: PerformanceWarning: Increasing number of chunks by factor of 37
  result = blockwise(

mitgcm = xr.merge([mitgcm, mixpods.pdf_N2S2(mitgcm)])
mitgcm = mixpods.load(mitgcm)
mitgcm

<xarray.Dataset>
Dimensions:                (depth: 100, time: 174000, RF: 101, N2T_bins: 29,
                            S2_bins: 29, enso_transition_phase: 7, stat: 2,
                            N2_bins: 29, Rig_T_bins: 9)
Coordinates: (12/20)
  * depth                  (depth) float32 -248.8 -246.2 -243.8 ... -3.75 -1.25
    latitude               float64 0.025
    longitude              float64 -140.0
  * time                   (time) datetime64[ns] 1998-12-31T18:00:00 ... 2018...
  * RF                     (RF) float64 -250.0 -247.5 -245.0 ... -5.0 -2.5 0.0
    eucmax                 (time) float32 dask.array<chunksize=(162,), meta=np.ndarray>
    ...                     ...
  * S2_bins                (S2_bins) object [-5.0, -4.9) ... [-2.200000000000...
  * enso_transition_phase  (enso_transition_phase) object 'none' ... 'all'
  * stat                   (stat) object 'mean' 'count'
  * N2_bins                (N2_bins) object [-5.0, -4.9) ... [-2.200000000000...
    bin_areas              (N2T_bins, S2_bins) float64 0.01 0.01 ... 0.01 0.01
  * Rig_T_bins             (Rig_T_bins) object (-1.6, -1.4000000000000001] .....
Data variables: (12/31)
    DFrI_TH                (depth, time) float32 dask.array<chunksize=(100, 6000), meta=np.ndarray>
    KPPdiffKzT             (depth, time) float32 dask.array<chunksize=(100, 6000), meta=np.ndarray>
    KPPg_TH                (depth, time) float32 dask.array<chunksize=(100, 6000), meta=np.ndarray>
    KPPhbl                 (time) float32 dask.array<chunksize=(6000,), meta=np.ndarray>
    KPPviscAz              (depth, time) float32 dask.array<chunksize=(100, 6000), meta=np.ndarray>
    SSH                    (time) float32 dask.array<chunksize=(6000,), meta=np.ndarray>
    ...                     ...
    Rig_T                  (RF, time) float32 dask.array<chunksize=(100, 6000), meta=np.ndarray>
    Rig                    (RF, time) float32 dask.array<chunksize=(100, 6000), meta=np.ndarray>
    eps                    (RF, time) float32 dask.array<chunksize=(100, 6000), meta=np.ndarray>
    n2s2pdf                (N2T_bins, S2_bins, enso_transition_phase) float64 ...
    eps_n2s2               (stat, N2_bins, S2_bins, enso_transition_phase) float64 ...
    eps_ri                 (stat, Rig_T_bins, enso_transition_phase) float64 ...
Attributes:
    easting:   longitude
    northing:  latitude
    title:     Station profile, index (i,j)=(1201,240)
mitgcm.u.cf.plot()
mitgcm.mldT.reset_coords(drop=True).cf.plot()
mitgcm.eucmax.reset_coords(drop=True).cf.plot()

[<matplotlib.lines.Line2D at 0x2ae9f91ebdc0>]
_images/mixpods-mom6_16_1.png 
mixpods.plot_n2s2pdf(mitgcm.n2s2pdf.sel(enso_transition_phase="none"))

<matplotlib.contour.QuadContourSet at 0x2ae9b69f4ac0>
_images/mixpods-mom6_17_1.png

Error#

shape = (101, 174000)
dims = ("RF", "time")
eps = mitgcm.eps.copy(
    data=dask.array.random.random_sample(shape, chunks=((100, 1), 8760))
)
N2 = mitgcm.N2T.copy(
    data=dask.array.random.random_sample(shape, chunks=((100, 1), 8760)), name="N2"
)
S2 = mitgcm.S2.copy(
    data=dask.array.random.random_sample(shape, chunks=((100, 1), 8760)), name="S2"
)
enso_transition = mitgcm.enso_transition

bins = np.arange(-5, -2.05, 0.1)
func = ["mean", "std"]

out = [
    flox.xarray.xarray_reduce(
        eps,
        N2,
        S2,
        enso_transition,
        func=f,
        expected_groups=(bins, bins, None),
        isbin=(True, True, False),
    )
    for f in func
]

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In [100], line 4
      2 dims = ("RF", "time")
      3 eps = mitgcm.eps.copy(data=dask.array.random.random_sample(shape, chunks=((100,1), 8760)))
----> 4 N2 = mitgcm.N2T.copy(data=dask.array.random.random_sample(shape, chunks=((100,1), 8760)), name="N2")
      5 S2 = mitgcm.S2.copy(data=dask.array.random.random_sample(shape, chunks=((100,1), 8760)))
      6 enso_transition = mitgcm.enso_transition

TypeError: DataArray.copy() got an unexpected keyword argument 'name'

rng = np.random.default_rng()

mitgcm.eps.shape

(101, 174000)

# shape = (10, 2000)
# chunks = ((shape[0]-1,1), 10)
shape = (101, 174000)
chunks = ((101,), 8760)
dims = ("RF", "time")
eps = xr.DataArray(
    data=dask.array.random.random_sample(shape, chunks=chunks), dims=dims
)
N2 = xr.DataArray(
    data=dask.array.random.random_sample(shape, chunks=chunks), dims=dims, name="N2"
)
S2 = xr.DataArray(
    data=dask.array.random.random_sample(shape, chunks=chunks), dims=dims, name="S2"
)
enso_transition = xr.DataArray(
    data=np.arange(5)[rng.integers(low=0, high=5, size=shape[-1])],
    dims=dims[1],
    name="enso",
)

bins = np.arange(-5, -2.05, 0.1)
func = ["mean", "std"]

out = [
    flox.xarray.xarray_reduce(
        eps,
        N2,
        S2,
        enso_transition,
        func=f,
        expected_groups=(bins, bins, None),
        isbin=(True, True, False),
    )
    for f in func
]

MOM6#

calculate ONI#

(
    oniobs.hvplot.line(x="time", label="obs")
    * onimom6.hvplot.line(x="time", label="MOM6")
).opts(ylabel="ONI [°C]")

oniobs.enso_transition_mask.plot()
onimom6.enso_transition_mask.plot(color="r")

[<matplotlib.lines.Line2D at 0x2b6d7fd20eb0>]
_images/mixpods-mom6_26_1.png

MOM6 sections#

Combine sections#
casename = "gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.001.mixpods"
mixpods.mom6_sections_to_zarr(casename)

/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask_jobqueue/core.py:255: FutureWarning: job_extra has been renamed to job_extra_directives. You are still using it (even if only set to []; please also check config files). If you did not set job_extra_directives yet, job_extra will be respected for now, but it will be removed in a future release. If you already set job_extra_directives, job_extra is ignored and you can remove it.
  warnings.warn(warn, FutureWarning)
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask_jobqueue/core.py:274: FutureWarning: env_extra has been renamed to job_script_prologue. You are still using it (even if only set to []; please also check config files). If you did not set job_script_prologue yet, env_extra will be respected for now, but it will be removed in a future release. If you already set job_script_prologue, env_extra is ignored and you can remove it.
  warnings.warn(warn, FutureWarning)
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask_jobqueue/pbs.py:82: FutureWarning: project has been renamed to account as this kwarg was used wit -A option. You are still using it (please also check config files). If you did not set account yet, project will be respected for now, but it will be removed in a future release. If you already set account, project is ignored and you can remove it.
  warnings.warn(warn, FutureWarning)
100%|██████████| 21/21 [02:28<00:00,  7.06s/it]
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/coding/times.py:360: FutureWarning: Index.ravel returning ndarray is deprecated; in a future version this will return a view on self.
  sample = dates.ravel()[0]
/glade/u/home/dcherian/pump/pump/mixpods.py:844: RuntimeWarning: Converting a CFTimeIndex with dates from a non-standard calendar, 'noleap', to a pandas.DatetimeIndex, which uses dates from the standard calendar.  This may lead to subtle errors in operations that depend on the length of time between dates.
  

casename = "gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.epbl.001.mixpods"
mixpods.mom6_sections_to_zarr(casename)

/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask_jobqueue/core.py:255: FutureWarning: job_extra has been renamed to job_extra_directives. You are still using it (even if only set to []; please also check config files). If you did not set job_extra_directives yet, job_extra will be respected for now, but it will be removed in a future release. If you already set job_extra_directives, job_extra is ignored and you can remove it.
  warnings.warn(warn, FutureWarning)
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask_jobqueue/core.py:274: FutureWarning: env_extra has been renamed to job_script_prologue. You are still using it (even if only set to []; please also check config files). If you did not set job_script_prologue yet, env_extra will be respected for now, but it will be removed in a future release. If you already set job_script_prologue, env_extra is ignored and you can remove it.
  warnings.warn(warn, FutureWarning)
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask_jobqueue/core.py:255: FutureWarning: job_extra has been renamed to job_extra_directives. You are still using it (even if only set to []; please also check config files). If you did not set job_extra_directives yet, job_extra will be respected for now, but it will be removed in a future release. If you already set job_extra_directives, job_extra is ignored and you can remove it.
  warnings.warn(warn, FutureWarning)
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask_jobqueue/core.py:274: FutureWarning: env_extra has been renamed to job_script_prologue. You are still using it (even if only set to []; please also check config files). If you did not set job_script_prologue yet, env_extra will be respected for now, but it will be removed in a future release. If you already set job_script_prologue, env_extra is ignored and you can remove it.
  warnings.warn(warn, FutureWarning)
100%|██████████| 31/31 [01:59<00:00,  3.85s/it]
/glade/u/home/dcherian/pump/pump/mixpods.py:846: RuntimeWarning: Converting a CFTimeIndex with dates from a non-standard calendar, 'noleap', to a pandas.DatetimeIndex, which uses dates from the standard calendar.  This may lead to subtle errors in operations that depend on the length of time between dates.
  mom6tao["time"] = mom6tao.indexes["time"].to_datetimeindex()

Read sections#

dask.config.set(scheduler=client)
m = mixpods.read_mom6_sections(
    "gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.001.mixpods"
)

/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask_jobqueue/core.py:255: FutureWarning: job_extra has been renamed to job_extra_directives. You are still using it (even if only set to []; please also check config files). If you did not set job_extra_directives yet, job_extra will be respected for now, but it will be removed in a future release. If you already set job_extra_directives, job_extra is ignored and you can remove it.
  warnings.warn(warn, FutureWarning)
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask_jobqueue/core.py:274: FutureWarning: env_extra has been renamed to job_script_prologue. You are still using it (even if only set to []; please also check config files). If you did not set job_script_prologue yet, env_extra will be respected for now, but it will be removed in a future release. If you already set job_script_prologue, env_extra is ignored and you can remove it.
  warnings.warn(warn, FutureWarning)
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask_jobqueue/pbs.py:82: FutureWarning: project has been renamed to account as this kwarg was used wit -A option. You are still using it (please also check config files). If you did not set account yet, project will be respected for now, but it will be removed in a future release. If you already set account, project is ignored and you can remove it.
  warnings.warn(warn, FutureWarning)
100%|██████████| 21/21 [02:14<00:00,  6.42s/it]
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/coding/times.py:360: FutureWarning: Index.ravel returning ndarray is deprecated; in a future version this will return a view on self.
  sample = dates.ravel()[0]
/glade/u/home/dcherian/pump/pump/mixpods.py:847: RuntimeWarning: Converting a CFTimeIndex with dates from a non-standard calendar, 'noleap', to a pandas.DatetimeIndex, which uses dates from the standard calendar.  This may lead to subtle errors in operations that depend on the length of time between dates.
  mom6tao["time"] = mom6tao.indexes["time"].to_datetimeindex()
/glade/u/home/dcherian/pump/pump/mixpods.py:859: UserWarning: Kv_v not present. Assuming equal to Kv_u
  warnings.warn("Kv_v not present. Assuming equal to Kv_u")

m.drop_vars(["average_DT", "average_T1", "average_T2", "time_bnds"]).chunk(
    {"time": 365 * 24}
).to_zarr(
    "/glade/scratch/dcherian/archive/gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.001.mixpods/ocn/moorings/tao.zarr",
    mode="w",
)

m.sel(yh=0, method="nearest")[["ePBL_h_ML", "mlotst"]].to_array().hvplot.line(
    by="variable", x="time"
)

m.sel(yh=0, method="nearest")[["Kd_heat", "Kd_ePBL"]].to_array().hvplot.line(
    by="variable", groupby="time", logy=True, ylim=(1e-6, 1e-1), xlim=(0, 500)
)

mom6140 = mixpods.load_mom6_sections(casename).load()

/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xgcm/grid.py:1471: UserWarning: Metric at ('time', 'zi', 'yh', 'xh') being interpolated from metrics at dimensions ('time', 'zl', 'yh', 'xh'). Boundary value set to 'extend'.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xgcm/grid.py:1471: UserWarning: Metric at ('time', 'zi', 'yh', 'xh') being interpolated from metrics at dimensions ('time', 'zl', 'yh', 'xh'). Boundary value set to 'extend'.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xgcm/grid.py:1471: UserWarning: Metric at ('time', 'zi', 'yq', 'xh') being interpolated from metrics at dimensions ('time', 'zl', 'yh', 'xh'). Boundary value set to 'extend'.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/cf_xarray/accessor.py:1638: UserWarning: Variables {'areacello'} not found in object but are referred to in the CF attributes.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xgcm/grid.py:1471: UserWarning: Metric at ('time', 'zi', 'yh', 'xh') being interpolated from metrics at dimensions ('time', 'zl', 'yh', 'xh'). Boundary value set to 'extend'.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xgcm/grid.py:1471: UserWarning: Metric at ('time', 'zi', 'yq', 'xh') being interpolated from metrics at dimensions ('time', 'zl', 'yh', 'xh'). Boundary value set to 'extend'.
  warnings.warn(

mom6140.uo.cf.plot(robust=True)
mom6140.eucmax.cf.plot()
mom6140.mldT.cf.plot(lw=0.5, color="k")

mixpods.plot_n2s2pdf(mom6140.n2s2pdf.sel(enso_transition_phase="none"))

<matplotlib.contour.QuadContourSet at 0x2acdb4c75e10>
_images/mixpods-mom6_39_1.png

Pick simulations#

mom6140 = mixpods.load_mom6_sections(
    "gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.001.mixpods"
)

/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xgcm/grid.py:1515: UserWarning: Metric at ('time', 'zi', 'yh', 'xh') being interpolated from metrics at dimensions ('time', 'zl', 'yh', 'xh'). Boundary value set to 'extend'.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xgcm/grid.py:1515: UserWarning: Metric at ('time', 'zi', 'yh', 'xh') being interpolated from metrics at dimensions ('time', 'zl', 'yh', 'xh'). Boundary value set to 'extend'.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xgcm/grid.py:1515: UserWarning: Metric at ('time', 'zi', 'yh', 'xh') being interpolated from metrics at dimensions ('time', 'zl', 'yh', 'xh'). Boundary value set to 'extend'.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xgcm/grid.py:1515: UserWarning: Metric at ('time', 'zi', 'yq', 'xh') being interpolated from metrics at dimensions ('time', 'zl', 'yh', 'xh'). Boundary value set to 'extend'.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/cf_xarray/accessor.py:1663: UserWarning: Variables {'areacello'} not found in object but are referred to in the CF attributes.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/coding/cftime_offsets.py:1130: FutureWarning: Argument `closed` is deprecated in favor of `inclusive`.
  return pd.date_range(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xgcm/grid.py:1515: UserWarning: Metric at ('time', 'zi', 'yh', 'xh') being interpolated from metrics at dimensions ('time', 'zl', 'yh', 'xh'). Boundary value set to 'extend'.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xgcm/grid.py:1515: UserWarning: Metric at ('time', 'zi', 'yq', 'xh') being interpolated from metrics at dimensions ('time', 'zl', 'yh', 'xh'). Boundary value set to 'extend'.
  warnings.warn(

epbl = mixpods.load_mom6_sections(
    "gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.epbl.001.mixpods"
)

/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xgcm/grid.py:1515: UserWarning: Metric at ('time', 'zi', 'yh', 'xh') being interpolated from metrics at dimensions ('time', 'zl', 'yh', 'xh'). Boundary value set to 'extend'.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xgcm/grid.py:1515: UserWarning: Metric at ('time', 'zi', 'yh', 'xh') being interpolated from metrics at dimensions ('time', 'zl', 'yh', 'xh'). Boundary value set to 'extend'.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xgcm/grid.py:1515: UserWarning: Metric at ('time', 'zi', 'yq', 'xh') being interpolated from metrics at dimensions ('time', 'zl', 'yh', 'xh'). Boundary value set to 'extend'.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/cf_xarray/accessor.py:1663: UserWarning: Variables {'areacello'} not found in object but are referred to in the CF attributes.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/coding/cftime_offsets.py:1130: FutureWarning: Argument `closed` is deprecated in favor of `inclusive`.
  return pd.date_range(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xgcm/grid.py:1515: UserWarning: Metric at ('time', 'zi', 'yh', 'xh') being interpolated from metrics at dimensions ('time', 'zl', 'yh', 'xh'). Boundary value set to 'extend'.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xgcm/grid.py:1515: UserWarning: Metric at ('time', 'zi', 'yq', 'xh') being interpolated from metrics at dimensions ('time', 'zl', 'yh', 'xh'). Boundary value set to 'extend'.
  warnings.warn(

datasets = {
    "TAO": tao_gridded,
    "MITgcm": mitgcm,
    "MOM6 KPP": mixpods.load(mom6140),
    "MOM6 ePBL": mixpods.load(epbl),
}
from datatree import DataTree

tree = DataTree.from_dict(
    {k: v.sel(time=slice("2000", "2017")) for k, v in datasets.items()}
)

tree["MOM6 KPP"].ds["time"] = tree["MOM6 KPP"].ds["time"] - pd.Timedelta("7h")

Verify depth is normalized#

for name, ds in datasets.items():
    (ds.cf["sea_water_x_velocity"].cf["Z"].plot(marker=".", ls="none", label=name))
plt.legend()

<matplotlib.legend.Legend at 0x2aed634274c0>
_images/mixpods-mom6_46_1.png

Compare EUC maximum and MLD#

Monthly climatology#

import tqdm

for node in tqdm.tqdm(tree.children):
    tree[node]["mldT"] = tree[node]["mldT"].load()

100%|██████████| 4/4 [00:52<00:00, 13.10s/it]

def to_dataset(tree):
    concat = xr.concat([v.ds for v in tree.children.values()], dim="node")
    concat["node"] = list(tree.children.keys())
    return concat


clim = to_dataset(
    tree.map_over_subtree(
        lambda x: x.reset_coords()[["mldT", "eucmax"]].groupby("time.month").mean()
    )
).load()
clim

<xarray.Dataset>
Dimensions:  (node: 4, month: 12)
Coordinates:
  * month    (month) int64 1 2 3 4 5 6 7 8 9 10 11 12
  * node     (node) <U9 'TAO' 'MITgcm' 'MOM6 KPP' 'MOM6 ePBL'
Data variables:
    mldT     (node, month) float64 -28.23 -18.01 -14.84 ... -25.31 -32.05 -33.92
    eucmax   (node, month) float64 -108.0 -104.2 -95.85 ... -119.1 -120.1 -120.9
hv.Layout(
    [
        hv.Overlay(
            [
                tree[node]
                .ds["mldT"]
                .reset_coords(drop=True)
                .groupby("time.month")[month]
                .hvplot.hist(label=node, legend=True)
                .opts(frame_width=150)
                for node in tree.children
            ]
        ).opts(title=str(month))
        for month in np.arange(1, 13)
    ]
).cols(4)

(clim.mldT.hvplot(by="node") + clim.eucmax.hvplot(by="node")).cols(1)

mixpods.plot_timeseries(tree, "mldT", obs="TAO")

datasets["TAO"].eucmax.attrs["long_name"] = "EUC maximum"
mixpods.plot_timeseries(tree.sel(time="2008"), "eucmax", obs="TAO")

MLD Maps#

def read_sfc(casename):
    path = f"/glade/scratch/dcherian/archive/{casename}/ocn/hist/*sfc*00[4-7]*.nc"
    sfc = xr.open_mfdataset(
        path, data_vars="minimal", coords="minimal", compat="override", parallel=True
    )
    return sfc


sfc = DataTree()
casenames = {
    "MOM6 KPP": "gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.001.mixpods",
    "MOM6 ePBL": "gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.epbl.001.mixpods",
}

for name, casename in casenames.items():
    sfc[name] = DataTree(read_sfc(casename))

sfc

<xarray.DatasetView>
Dimensions:  ()
Data variables:
    *empty*
mldclim = sfc.map_over_subtree(
    lambda node: node.mlotst.groupby("time.month").mean()
).compute()

filepath = (
    "/glade/work/gmarques/cesm/datasets/MLD/deBoyer/deBoyer_MLD_remapped_to_tx06v1.nc"
)
mldclim["obs"] = DataTree(xr.open_dataset(filepath))
mldclim

/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/coding/times.py:699: SerializationWarning: Unable to decode time axis into full numpy.datetime64 objects, continuing using cftime.datetime objects instead, reason: dates out of range
  dtype = _decode_cf_datetime_dtype(data, units, calendar, self.use_cftime)
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/core/indexing.py:524: SerializationWarning: Unable to decode time axis into full numpy.datetime64 objects, continuing using cftime.datetime objects instead, reason: dates out of range
  return np.asarray(array[self.key], dtype=None)

<xarray.DatasetView>
Dimensions:  ()
Data variables:
    *empty*
mldclim.

mldclim["MOM6 KPP"] = DataTree(
    read_sfc("gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.001.mixpods")
    .mlotst.groupby("time.month")
    .mean()
)

mldclim["MOM6 ePBL"] = DataTree(
    read_sfc("gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.epbl.001.mixpods")
    .mlotst.groupby("time.month")
    .mean()
)

(
    (mlotst - mldclim.mld.groupby("time.month").first())
    .sel(yh=slice(-5, 5), xh=slice(-170, -80))
    .plot(col="month", robust=True, col_wrap=3, aspect=2)
)

<xarray.plot.facetgrid.FacetGrid at 0x2b2239f16c50>
_images/mixpods-mom6_62_1.png

Mixing layer depth?#

todrop = ["TAO"]
notao = DataTree.from_dict(
    {k: v.ds for k, v in tree.children.items() if k not in todrop}
)
with cfxr.set_options(custom_criteria=mixing_layer_criteria):
    plot = mixpods.plot_timeseries(notao, "boundary_layer_depth", obs=None)
plot

/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/cf_xarray/accessor.py:1663: UserWarning: Variables {'areacello'} not found in object but are referred to in the CF attributes.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/cf_xarray/accessor.py:1663: UserWarning: Variables {'areacello'} not found in object but are referred to in the CF attributes.
  warnings.warn(

Mean profiles: mean, std#

Model mean S2 is lower by a factor of 2

tree = tree.assign(
    dTdz=lambda ds: ds.cf["sea_water_potential_temperature"]
    .reset_coords(drop=True)
    .cf.differentiate("Z")
    .load()
)
tree = tree.assign(
    dSdz=lambda ds: ds.cf["sea_water_salinity"]
    .reset_coords(drop=True)
    .cf.differentiate("Z")
    .load()
)
tree = tree.assign({"4N2T": lambda ds: 4 * ds.N2T.load()})

/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/cf_xarray/accessor.py:1663: UserWarning: Variables {'areacello'} not found in object but are referred to in the CF attributes.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/cf_xarray/accessor.py:1663: UserWarning: Variables {'areacello'} not found in object but are referred to in the CF attributes.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/cf_xarray/accessor.py:1663: UserWarning: Variables {'areacello'} not found in object but are referred to in the CF attributes.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/cf_xarray/accessor.py:1663: UserWarning: Variables {'areacello'} not found in object but are referred to in the CF attributes.
  warnings.warn(

kwargs = dict(
    show_grid=True,
    title="",
    invert_axes=True,
    legend_position="right",
    frame_width=150,
    frame_height=500,
)

limkwargs = dict(
    ylim=(-1e-3, 1.9e-3),
    xlim=(-250, 0),
)


def plot_profile_fill(datasets, var, label):
    from datatree import DataTree

    if isinstance(datasets, DataTree):
        datasets = datasets.children
    return mixpods.map_hvplot(
        lambda ds, name: mixpods.plot_profile_fill(ds.ds.cf[var].load(), name), datasets
    ).opts(**kwargs, ylabel=label)


S2 = plot_profile_fill(tree, "S2", "S^2")
N2 = plot_profile_fill(tree, "N2T", "N_T^2")
Ri = (
    (
        mixpods.map_hvplot(
            lambda ds, name: mixpods.cfplot(ds.Rig_T.load().median("time"), name),
            datasets,
        )
        * hv.HLine(0.25).opts(color="k", line_width=0.5)
    )
    .opts(**kwargs, **limkwargs, ylabel="median Ri_g^T")
    .opts(ylim=(0, 1))
)

(S2 + N2 + Ri).opts(sizing_mode="stretch_width")

T = plot_profile_fill(tree, "sea_water_potential_temperature", "T")
S = plot_profile_fill(tree, "sea_water_salinity", "S")
u = plot_profile_fill(tree, "sea_water_x_velocity", "u")
# dTdz = plot_profile_fill(tree, "dTdz", "∂T/∂z")
# dSdz = plot_profile_fill(tree, "dSdz", "∂S/∂z")

/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/cf_xarray/accessor.py:1663: UserWarning: Variables {'areacello'} not found in object but are referred to in the CF attributes.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/cf_xarray/accessor.py:1663: UserWarning: Variables {'areacello'} not found in object but are referred to in the CF attributes.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/cf_xarray/accessor.py:1663: UserWarning: Variables {'areacello'} not found in object but are referred to in the CF attributes.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/cf_xarray/accessor.py:1663: UserWarning: Variables {'areacello'} not found in object but are referred to in the CF attributes.
  warnings.warn(

(T + S + u.opts(ylim=(-0.5, 1.5))).opts(sizing_mode="stretch_width")

Remap to EUC coordinate#

gcmeq.coords["zeuc"] = gcmeq.depth - gcmeq.eucmax
remapped = flox.xarray.xarray_reduce(
    gcmeq.drop_vars(["SSH", "KPPhbl", "mld", "eucmax"], errors="ignore"),
    "zeuc",
    dim="depth",
    func="mean",
    expected_groups=(np.arange(-250, 250.1, 5),),
    isbin=(True,),
)
remapped["zeuc_bins"].attrs["units"] = "m"
remapped

<xarray.Dataset>
Dimensions:          (time: 174000, longitude: 4, zeuc_bins: 100)
Coordinates:
    latitude         float64 0.025
  * longitude        (longitude) float64 -155.0 -140.0 -125.0 -110.0
  * time             (time) datetime64[ns] 1998-12-31T18:00:00 ... 2018-11-06...
    cool_mask        (time) bool True True True True ... False False False False
    warm_mask        (time) bool False False False False ... True True True True
  * zeuc_bins        (zeuc_bins) object (-250.0, -245.0] ... (245.0, 250.0]
Data variables: (12/25)
    DFrI_TH          (time, longitude, zeuc_bins) float32 dask.array<chunksize=(6000, 1, 100), meta=np.ndarray>
    KPPdiffKzT       (time, longitude, zeuc_bins) float32 dask.array<chunksize=(6000, 1, 100), meta=np.ndarray>
    KPPg_TH          (time, longitude, zeuc_bins) float32 dask.array<chunksize=(6000, 1, 100), meta=np.ndarray>
    KPPviscAz        (time, longitude, zeuc_bins) float32 dask.array<chunksize=(6000, 1, 100), meta=np.ndarray>
    VISrI_Um         (time, longitude, zeuc_bins) float32 dask.array<chunksize=(6000, 1, 100), meta=np.ndarray>
    VISrI_Vm         (time, longitude, zeuc_bins) float32 dask.array<chunksize=(6000, 1, 100), meta=np.ndarray>
    ...               ...
    S2               (time, longitude, zeuc_bins) float32 dask.array<chunksize=(6000, 1, 100), meta=np.ndarray>
    shear            (time, longitude, zeuc_bins) float32 dask.array<chunksize=(6000, 1, 100), meta=np.ndarray>
    N2               (time, longitude, zeuc_bins) float32 dask.array<chunksize=(6000, 1, 100), meta=np.ndarray>
    shred2           (time, longitude, zeuc_bins) float32 dask.array<chunksize=(6000, 1, 100), meta=np.ndarray>
    Ri               (time, longitude, zeuc_bins) float32 dask.array<chunksize=(6000, 1, 100), meta=np.ndarray>
    enso_transition  (zeuc_bins, time) <U12 'La-Nina cool' ... 'El-Nino warm'
Attributes:
    easting:   longitude
    northing:  latitude
    title:     Station profile, index (i,j)=(1201,240)
remapped = remapped.persist()

cluster.scale(6)

Seasonal median Ri: (0, 140)#

remapped["longitude"] = [-155.0, -140, -125, -110]

fg = (
    remapped.Ri.groupby("time.season")
    .median()
    .sel(season=["DJF", "MAM", "JJA", "SON"], longitude=[-140, -110])
    .cf.plot(
        col="season", row="longitude", xlim=(0, 1.5), ylim=(-20, None), label="MITgcm"
    )
)
fg.data = tao_Ri.cf.sel(quantile=0.5, longitude=[-140, -110])
fg.map_dataarray_line(
    xr.plot.line, x=None, y="zeuc", hue="season", color="k", label="TAO"
)
fg.map(lambda: plt.axvline(0.25))
fg.axes[-1, -1].legend(bbox_to_anchor=(1.5, 1))

<matplotlib.legend.Legend at 0x2add93b58220>
_images/mixpods-mom6_78_1.png

Stability diagram: 4N2-S2 plane#

Warner & Moum (2019):

Both velocity and temperature are hourly averaged before interpolation to 5-m vertical bins

Contours enclose 50% of data

mixpods.plot_stability_diagram_by_dataset(datasets)
_images/mixpods-mom6_80_0.png

Questions:

  1. internal wave spectrum?

  2. minor axis is smaller for model, particularly MITgcm-> too diffusive?

  3. Can we think of this as the model relaxing to critical Ri too quickly

mixpods.plot_stability_diagram_by_phase(datasets, obs="TAO", fig=None)
_images/mixpods-mom6_82_0.png

Composed#

excluding MLD#

  • have not matched vertical grid spacing yet.

  • minor axis is smaller for increasing shear with KPP

  • prandtl number

  • interpolate to TAO grid

  • hybrid HYCOM-like case may have equatorial refinement of vertical grid above EUC

  • compare hybrid coordinate case resolution

  • instrumental error on vertical gradients for S2, N2;

    • if added to model, would it change the bottom tail.

  • Horizontal or vertical res:

    • MOM6 at 5km? hah

    • TIW compositing? La-NIna dec?

      • top end’s are similar, so maybe not?

  • sensitivity of cold tongue mean state

    • horizontal viscosity; controls EUC strength? (increase?)

    • vertical res; yan li’s paper

    • POP :thermocline always too diffuse and EUC too thick;

    • higher order interpolants can control amount diffusion when remapping

  • get S2, N2 from TPOSE?

fig = plt.figure(constrained_layout=True, figsize=(12, 6))

subfigs = fig.subfigures(2, 1, height_ratios=[1.2, 1])
top = subfigs[0].subfigures(1, 3, width_ratios=[1, 10, 1])

mixpods.plot_stability_diagram_by_dataset(datasets, fig=top[1])
mixpods.plot_stability_diagram_by_phase(datasets, fig=subfigs[1])
_images/mixpods-mom6_84_0.png

including MLD#

fig = plt.figure(constrained_layout=True, figsize=(12, 6))
subfigs = fig.subfigures(2, 1, height_ratios=[1.2, 1])
top = subfigs[0].subfigures(1, 3, width_ratios=[1, 5, 1])

mixpods.plot_stability_diagram_by_dataset(datasets, fig=top[1])
mixpods.plot_stability_diagram_by_phase(datasets, fig=subfigs[1])
_images/mixpods-mom6_86_0.png 
(
    tao_gridded.n2s2pdf.sel(
        enso_transition_phase=[
            "La-Nina cool",
            "La-Nina warm",
            "El-Nino warm",
            "El-Nino cool",
        ]
    )
    .sum("N2_bins")
    .plot(hue="enso_transition_phase")
)
plt.figure()
(
    tao_gridded.n2s2pdf.sel(
        enso_transition_phase=[
            "La-Nina cool",
            "La-Nina warm",
            "El-Nino warm",
            "El-Nino cool",
        ]
    )
    .sum("S2_bins")
    .plot(hue="enso_transition_phase")
)

[<matplotlib.lines.Line2D at 0x2b36e2bc2ad0>,
 <matplotlib.lines.Line2D at 0x2b36e2bc32e0>,
 <matplotlib.lines.Line2D at 0x2b36e2bc3310>,
 <matplotlib.lines.Line2D at 0x2b36e0ca52d0>]
_images/mixpods-mom6_87_1.png _images/mixpods-mom6_87_2.png 
oni = pump.obs.process_oni().sel(time=slice("2005-Oct", "2015"))
(
    oni.hvplot.step(color="k")
    + pump.obs.make_enso_transition_mask()
    .sel(time=slice("2005-Jun", "2015"))
    .reset_coords(drop=True)
    .hvplot.step()
).cols(1)

ε Ri TAO check#

data = tao_gridded.drop_vars("zeuc")
epsZ = data.eps.cf["Z"].cf.sel(Z=slice(-69, -29))

# χpod data are available at a subset of depths
Ri = np.log10(mixpods.reindex_Z_to(data.Rig_T, epsZ))
# Ri_bins = np.logspace(np.log10(0.025), np.log10(2), 11)
Ri_bins = np.arange(-1.6, 0.4, 0.2)  # in log space from Sally
# Ri_bins = np.array([0.02, 0.04, 0.06, 0.1, 0.3, 0.5, 0.7, 1.5, 2])
Ri_kwargs = {
    "expected_groups": (Ri_bins, None),
    "isbin": (True, False),
}

mpl.rcParams["figure.dpi"] = 120

(
    data.eps.count("depthchi").hvplot.line(x="time")
    + Ri.count("depthchi").hvplot.line(x="time")
).cols(1)

data = tao_gridded.eps_n2s2.sel(
    enso_transition_phase=[
        "El-Nino warm",
        "La-Nina warm",
        "El-Nino cool",
        "La-Nina cool",
    ]
)

fg = data.to_dataset("stat").plot.scatter(
    x="S2_bins",
    y="N2_bins",
    hue="mean",  # markersize="count",
    col="enso_transition_phase",
    col_wrap=2,
    norm=mpl.colors.LogNorm(1e-7, 1e-6),
    cmap=mpl.cm.Greys,
    robust=True,
    xlim=(-4.5, -2.5),
    ylim=(-4.5, -2.5),
)
fg.set_titles("{value}")

np.log10(tao_gridded.eps_ri).sel(stat="mean").plot(
    hue="enso_transition_phase", ylim=(-7.5, -6)
)

/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/plot/plot.py:454: UserWarning: The handle <matplotlib.lines.Line2D object at 0x2aea35fc2b00> has a label of '____________' which cannot be automatically added to the legend.
  ax.legend(handles=primitive, labels=list(hueplt.to_numpy()), title=hue_label)

[<matplotlib.lines.Line2D at 0x2aea4088c520>,
 <matplotlib.lines.Line2D at 0x2aea4088fd90>,
 <matplotlib.lines.Line2D at 0x2aea4088ee90>,
 <matplotlib.lines.Line2D at 0x2aea4088e6e0>,
 <matplotlib.lines.Line2D at 0x2aea4088eec0>,
 <matplotlib.lines.Line2D at 0x2aea35fc2b00>,
 <matplotlib.lines.Line2D at 0x2aea35fc39a0>]
_images/mixpods-mom6_94_2.png

Turbulence#

Histograms: Ri, ε#

for ds in datasets.values():
    ds.Rig_T.load()

from xhistogram.xarray import histogram

bins = np.linspace(-11, -4, 101)
handles = [
    histogram(
        np.log10(subtree.ds.eps.cf.sel(Z=slice(-69, -29))), bins=bins, density=True
    ).hvplot.step(label=name)
    for name, subtree in tree.children.items()
]

hv.Overlay(
    handles
)  # .opts(opts.Histogram(alpha=0.5, fill_alpha=0, hover_fill_alpha=1))

bins = np.linspace(-0.5, 1.5, 61)
handles = [
    histogram(
        ds.Rig_T.reset_coords(drop=True).cf.sel(Z=slice(-69, -29)),
        density=True,
        bins=bins,
    ).hvplot.step(label=name)
    for name, ds in datasets.items()
]
(hv.Overlay(handles) * hv.VLine(0.25).opts(line_color="k")).opts(ylabel="PDF")

ε vs Ri#

Possible reasons for difference

  1. My MLD seems deeper even though I use the same ΔT=0.15 threshold. It could be that they’ve used Akima splines to interpolate in the vertical

  2. ~I’ve filtered to DCL, so accounting for MLD and EUC movement. I’m not sure they did that.~

Only 𝜒pods between 29 and 69 m are used in this analysis as deeper 𝜒pods are more strongly influenced by the variability of zEUC than by surface forcing.

TODO

  • Add \(ε_χ\) for MOM6

  • EUC strength is proportional to horizontal visc for 1° models

  • Do for K

  • composite by TIW strength

  • make composites off the equator: look at strong off-equatoirial du/dx; du/dz

f, ax = plt.subplots(
    1, len(datasets), sharex=True, sharey=True, constrained_layout=True
)
for axx, (name, ds) in zip(ax, datasets.items()):
    da = np.log10(ds.eps_ri)
    # da["Rig_T_bins"] = pd.Index(
    #    [pd.Interval(np.log10(i.left), np.log10(i.right))
    #     for i in tao_gridded.Rig_T_bins.data]
    # )
    (
        da.sel(stat="mean")
        .sel(enso_transition_phase=["La-Nina cool", "El-Nino warm"])
        .plot(hue="enso_transition_phase", marker=".", ax=axx, add_legend=False)
    )
    axx.set_title(name)

    ticks = [0.04, 0.1, 0.25, 0.63, 1.6]
    # axx.set_yticks([-7.25, -7, -6.5, -6])
    axx.set_ylim([-7.5, -5.5])
    axx.set_xticks(np.log10(ticks))
    axx.set_xticklabels([f"{a:.2f}" for a in ticks])
    # axx.tick_params(axis='x', rotation=30)
ax[0].set_ylabel("ε")
dcpy.plots.linex(np.log10(0.25), ax=ax)
dcpy.plots.clean_axes(ax)
f.set_size_inches((8, 2))
_images/mixpods-mom6_103_0.png 
np.log10(tao_gridded.eps_n2s2).plot(robust=True, col="enso_transition_phase")

<xarray.plot.facetgrid.FacetGrid at 0x2acdb3799ae0>
_images/mixpods-mom6_104_1.png

Seasonal mean heat flux#

remapped.Jq.load()

<xarray.DataArray 'Jq' (time: 174000, zeuc: 100)>
array([[        nan,         nan, -0.07841657, ...,         nan,
                nan,         nan],
       [        nan,         nan, -0.07973828, ...,         nan,
                nan,         nan],
       [        nan,         nan, -0.08094554, ...,         nan,
                nan,         nan],
       ...,
       [        nan, -0.07447129,         nan, ...,         nan,
                nan,         nan],
       [        nan, -0.07471326,         nan, ...,         nan,
                nan,         nan],
       [        nan, -0.07509062,         nan, ...,         nan,
                nan,         nan]])
Coordinates:
    latitude   float64 0.025
    longitude  float64 -140.0
  * time       (time) datetime64[ns] 1998-12-31T18:00:00 ... 2018-11-06T17:00:00
  * zeuc       (zeuc) float64 -247.5 -242.5 -237.5 -232.5 ... 237.5 242.5 247.5
import hvplot.xarray

(
    remapped.Jq.groupby("time.season")
    .mean()
    .reindex(season=["DJF", "MAM", "JJA", "SON"])
    .cf.plot.line(col="season")
)

<xarray.plot.facetgrid.FacetGrid at 0x2affc38316c0>
_images/mixpods-mom6_108_1.png

Test out available variables#

  1. KPP_Kheat is non-zero only in KPP_OBL

(mom6140.Tflx_dia_diff * 1025 * 4000).cf.plot(robust=True)

<matplotlib.collections.QuadMesh at 0x2ba58a988fa0>
_images/mixpods-mom6_110_1.png 
mom6140.Kv_u.cf.plot(norm=mpl.colors.LogNorm(1e-6, 1e-2))

<matplotlib.collections.QuadMesh at 0x2ba58a007cd0>
_images/mixpods-mom6_111_1.png 
mom6140.h.median("time").plot()
mom6140.h.mean("time").plot()

[<matplotlib.lines.Line2D at 0x2ba58a547c10>]
_images/mixpods-mom6_112_1.png 
mom6140.Kd_heat.cf.plot(norm=mpl.colors.LogNorm(1e-6, 1e-2))
mom6140.KPP_OBLdepth.plot(color="r")

plt.figure()
mom6140.KPP_kheat.cf.plot(norm=mpl.colors.LogNorm(1e-6, 1e-2))
mom6140.KPP_OBLdepth.plot(color="r")

plt.figure()
(mom6140.KPP_kheat - mom6140.Kd_heat).cf.plot(
    robust=True, cmap=mpl.cm.mpl.cm.Reds_r, vmax=0
)

<matplotlib.collections.QuadMesh at 0x2ba58a8c38b0>
_images/mixpods-mom6_113_1.png _images/mixpods-mom6_113_2.png _images/mixpods-mom6_113_3.png

Experiment with manual combining#

import intake
from intake.source.utils import reverse_format

years = []
tiles = []
for file in files[:10]:
    p = pathlib.Path(file)
    params = reverse_format("__{year}.nc.{tile}", p.stem + p.suffix)
    years.append(params["year"])
    tiles.append(params["tile"])
years, tiles

import toolz as tlz

bytile = {}
for tile, v in tlz.groupby(tileid, files).items():
    bytile[tile] = xr.concat(read_raw_files(v, parallel=True), dim="time")


print("\n".join([hash(ds.yh.data.tolist()) for ds in list(bytile.values())]))

from functools import partial


def hash_coords(ds, axis):
    dims = ds.cf.axes[axis]
    data = np.concatenate([ds[dim].data for dim in dims])
    return hash(tuple(data))


grouped = bytile

for axis, concat_axis in [("X", "Y"), ("Y", "X")]:
    grouped = tlz.groupby(partial(hash_coords, axis=axis), grouped.values())
    grouped = {
        k: cfconcat(v, axis=concat_axis, join="exact") for k, v in grouped.items()
    }
    break
grouped

cfconcat(list(grouped.values()), "X")

combined = xr.combine_by_coords(list(bytile.values()), combine_attrs="override")


def raise_if_bad_index(combined):
    bad = []
    for idx in combined.indexes:
        index = combined.indexes[idx]
        if not index.is_monotonic or index.has_duplicates:
            bad.append(idx)
    if bad:
        raise ValueError(f"Indexes {idx} are either not monotonic or have duplicates")


def tileid(path):
    p = pathlib.Path(path)
    # print(p.suffix)
    # params = reverse_format("__{year}.nc.{tile}", p.stem + p.suffix)
    return int(p.suffix[1:])  # params["tile"]
    # years.append(params["year"])
    # tiles.append(params["tile"])


sorted_files = sorted(files, key=tileid)
for tile, files in groupby(sorted_files, tileid):
    print(tile, len(list(files)))

Test out ONI#

Phase labeling#

oni = pump.obs.process_oni().sel(time=slice("2005-Sep", None))
enso_transition = mixpods.make_enso_transition_mask(oni)
mixpods.plot_enso_transition(oni, enso_transition)

Replicate ONI calculation#

close enough!

ersst = xr.tutorial.open_dataset("ersstv5")

monthlyersst = (
    ersst.sst.cf.sortby("Y")
    .cf.sel(latitude=slice(-5, 5), longitude=slice(360 - 170, 360 - 120))
    .cf.mean(["X", "Y"])
    .resample(time="M")
    .mean()
    .load()
)

expected = mixpods.calc_oni(monthlyersst)

actual = pump.obs.process_oni()
actual.plot()
expected.plot()

(actual - expected).plot()

[<matplotlib.lines.Line2D at 0x2b6d2ec04580>]
_images/mixpods-mom6_119_1.png

Comparing to Warner and Moum code#

chipod.eps.sel(time="2007-01-14 10:35:00", method="nearest").load()

<xarray.DataArray 'eps' (depth: 5)>
nan 3.929e-07 2.846e-07 2.234e-07 4.491e-07
Coordinates:
    time        datetime64[ns] 2007-01-14T11:00:00
  * depth       (depth) float64 -69.0 -59.0 -49.0 -39.0 -29.0
    timeSeries  (depth) float64 69.0 59.0 49.0 39.0 29.0
    lat         (depth) float64 0.0 0.0 0.0 0.0 0.0
    lon         (depth) float64 -140.0 -140.0 -140.0 -140.0 -140.0
Attributes:
    long_name:              turbulence dissipation rate
    standard_name:          specific_turbulent_kinetic_energy_dissipation_in_...
    ncei_name:              turbulent kinetic energy (TKE) dissipation rate
    units:                  W kg-1
    FillValue:              -9999
    valid_min:              1e-12
    valid_max:              9.999999999999999e-06
    coverage_content_type:  physicalMeasurement
    grid_mapping:           crs
    source:                 inferred from fast thermistor spectral scaling
    references:             Moum J.N. and J.D. Nash, Mixing measurements on a...
    cell_methods:           depth: point, time: mean
    platform:               mooring
    instrument:             chipod
tao_gridded.enso_transition.loc[{"time": slice("2015-12", "2016-01")}]

<xarray.DataArray 'enso_transition' (time: 1488)>
'El-Nino warm' 'El-Nino warm' 'El-Nino warm' ... 'El-Nino warm' 'El-Nino warm'
Coordinates: (12/13)
    deepest             (time) float64 -300.0 -300.0 -300.0 ... -300.0 -300.0
    eucmax              (time) float64 -135.0 -140.0 -135.0 ... -140.0 -135.0
    latitude            float32 0.0
    longitude           float32 -140.0
    mld                 (time) float64 -9.0 -8.0 -8.0 -8.0 ... nan nan nan nan
    mldT                (time) float64 -10.0 -8.0 -8.0 -8.0 ... nan nan nan nan
    ...                  ...
    shallowest          (time) float64 -5.0 -5.0 -5.0 -5.0 ... -10.0 -10.0 -10.0
  * time                (time) datetime64[ns] 2015-12-01 ... 2016-01-31T23:00:00
    oni                 (time) float32 2.53 2.53 2.53 2.53 ... 2.53 2.53 2.53
    warm_mask           (time) bool True True True True ... False False False
    cool_mask           (time) bool False False False False ... True True True
    enso_transition     (time) <U12 'El-Nino warm' ... 'El-Nino warm'
Attributes:
    description:  Warner & Moum (2019) ENSO transition phase; El-Nino = ONI >...

Extra code#